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I.  INTRODUCTION 

For  many  years  it  has  been  observed  that  liquid-filled  projectiles  have  a 
proclivity  for  unusual  flight  behavior,  often  being  unstable  even  though  the 
same  projectile  with  a  solid  payload  is  stable.  The  history  of  various 
aspects  of  this  area  of  ballistics,  summaries  and  critiques  of  previous 

theoretical  work  and  a  discussion  of  experimental  techniques  are  given  in 
References  1  and  2.  The  distinction  between  the  steady  state  and  spin-up 
problems  is  basic. Only  the  steady  state  problem  for  an  incompressible 
fluid  in  a  cylinder  is  considered  here;  therefore,  the  undisturbed  basic  flow 
is  solid  body  rotation.  In  principle,  the  same  techniques  could  be  used  for 
the  spin-up  problem. 

To  put  this  work  in  perspective,  it  is  helpful  to  distinguish  between 

analytical  and  finite  difference  approaches  (see  References  1  and  2).  So  far, 
analytical  approaches  have  been  confined  to  linearized  problems.  Most  of  these 
are  asymptotic  for  Re  *  »;  the  case  Re  -*•  0  was  considered  by  Sedney  (unpub¬ 
lished,  see  Reference  1  for  a  summary).  The  Reynolds  number.  Re,  is  defined 
below.  Herbert^  analyzed  experimental  data  at  "low  Re"  and  constructed  a 
model  for  an  infinite  cylinder  from  which  he  deduced  results  for  1  <  Re  <  10'*. 

His  results  for  roll  moment  agree  with  experiment.  In  principle,  the  analyti¬ 
cal  approach  used  in  this  work  is  valid  for  all  Re;  it  is  "exact"  in  the  sense 

that  the  accuracy  of  the  results  can  be  improved  by  adding  more  eigenfunctions 
in  the  expansion.  In  practice,  treating  large  Re,  say  Re  >  1,000,  becomes 
tedious  using  the  techniques  presented  here. 

An  outline  of  the  process  in  going  from  the  Navier-Stokes  equations  to 
the  linear  partial  differential  equations  and  the  boundary  conditions  is  given 
now;  the  details  are  in  Reference  4. 

1.  The  Navier-Stokes  equations  are  written  in  an  inertial  reference  frame 
using  cylindrical  polar  coordinates  (r,e,x)  with  corresponding  velocity 
components  (u,v,w);  an  inertial  frame  is  used  even  though  there  is  an 
advantage  to  using  an  aerobal 1 istic  reference  frame  in  that  the  boundary 
conditions  are  easier  to  apply  (see  Reference  1). 

2.  To  obtain  boundary  conditions,  the  motion  of  the  cylinder  walls  is  deter¬ 
mined  from  the  projectile  motion  which  is  proportional  to  exp  i(ft-e)  where  t 
is  time  and  f  is,  in  general,  the  non-dimensional  complex  frequency  of  the 
projectile  motion;  for  convenience,  f  is  taken  to  be  real  here,  so  that  f  =  t, 
the  non-dimensional  coning  frequency  for  pure  coning  motion. 

3.  The  angle  of  the  coning  motion,  Kq,  is  assumed  small  and  only  linear 

departures  from  solid  body  rotation  are  considered.  The  velocity  components 
are  then  (-Kq  u*,  r  -Kq  v*,  -Kq  w*)  and  pressure  is  1/2  r2  -  Kq  P*.  The  * 

denotes  perturbation. 

4.  The  motion  of  the  cylinder  is  specified  and  the  response  of  the  fluid 
determined;  no  feedback  from  the  fluid  to  the  projectile  is  allowed. 

5.  It  is  convenient  to  introduce  complex  velocities  and  pressure,  denoted  by 
sub  _  and  defined  by 


cscscwwww 


(u*,v*,w*,p*)  =  Real  I  (u  ,v  ,w,p)  exp  i(ft-fne)| 


(1.1) 


where  the  azimuthal  wave  number,  m,  is  an  integer  (*)  but  only  m  =  1  is  needed 
for  the  linear  forced  motion  problem. 

The  resulting  linearized,  non-dimensional  Navier- Stokes  equations  are: 


i(f-m)  u-2v  =  -p  +  Re"^  [v^u  -  — ^  ^  ^  u  +  1— lli  y] 

■  ■  -  ■  r2  •  r2  ■ 


i  (f-m)  V  +  2u  =  (im/r)  p  +  Re'^  [v2v  -  -  —  v  -  -  u] 

r2  -  r2  - 


1  r„o..  m2 


i  (f-m)  w  =  -p  +  Re’  [72w - w] 


(ru)^  -  im  V  +  rw^  =  0  , 


(1.2a,b,c,d) 


where  72  =  3^2  +  (1/r)  3„  +  3  2  is  the  Laplacian  and  subscripts  denote 
partial  derivatives.  The  non-dimensional  variables  are  defined  after  (2.1), 

From  (1.2)  and  appropriate  boundary  conditions,  the  eigenvalue  and  forced 
motion  problems  are  defined  in  Section  II,  Particular  solutions  of  the  par¬ 
tial  differential  equations  which  simplify  application  of  boundary  conditions 
are  also  introduced  there,  and  finally,  the  eigenfunction  expansion  of  the 
solution.  The  latter  is  the  basis  of  the  spatial  eigenvalue  method  which  was 
used  by  B1  ennerhasset  and  Hall  5  in  their  study  of  the  Taylor  vortex  problem 
for  finite  length  cylinders. 

Methods  for  solving  the  eigenvalue  problem  are  discussed  in  Section 
III.  Obtaining  the  eigenvalues  is  central  to  the  method  and  is  a  nontrivial 
problem.  The  complex  eigenvalues  are  defined  by  a  6th  order  complex  system  of 
ordinary  differential  equations  which  is  not  self-adjoint  and  are  computed  by 
an  iterative  process.  Sufficiently  accurate  initial  estimates,  or  first 
guesses,  are  required  for  the  iteration  process  to  converge;  these  are  ob¬ 
tained  using  asymptotic  approximations  for  Re  -*■  0  and  Re  ♦  »,  variation  on  the 
axi Symmetric  Taylor  problem  for  Re  not  too  large,  or  extrapolation.  Addition¬ 
ally,  an  algoritln  was  developed  to  compute  an  initial  estimate  from  a 
previous  eigenvalue  by  a  trivial  arithmetic  operation.  Once  the  eigenvalues 
and  eigenfunctions  are  known,  the  coefficients  in  the  expansion  must  be  com¬ 
puted;  collocation  and  least  squares  methods  are  used. 

Formulas  for  pressure  and  moment  coefficients  are  presented  and  discussed 
in  Section  IV.  Asymptotic  forms  for  the  pressure  coefficient  for  f  +  0  and 
f  -►  00  are  al  so  given. 


Results  are  presented  in  Section  V.  Although  the  velocities  and  pressure 
are  available  in  the  output,  they  are  usually  not  extracted.  Instead,  the 
variation  of  pressure  and  moment  coefficients  with  Re,  f  or  aspect  ratio  are 
shown  since  these  are  the  items  of  primary  interest.  Some  comparisons  are 
made  with  results  from  experiments.  Several  comparisons  are  also  made  with 
results  from  the  finite  difference  method  of  Stri kwerda ,6  as  computed  by 
Nusca. Rather  limited  comparisons  are  made  with  results  from  other  numerical 
methods. 

Because  of  the  emphasis  given  results  from  the  Strikwerda^  method  in 
Section  V,  it  is  apropos  to  consider  some  features  of  it.  Since  it  solves  the 
nonlinear  Navier-Stokes  equations,  it  can  handle  large  K^;  and  it  computes 

both  the  roll  and  side  moment  coefficients.  At  present,  it  has  been  applied 
to  Re  <  100,  Reference  7.  CPU  time  on  the  VAX  8600  is  of  the  order  of  hours, 

increasing  with  Re.  For  example,  if  Re  =  10,  f  =  0.05,  and  aspect  ratio 
(length/diameter)  of  the  cylinder  is  3.0,  it  is  1.1  hours. 

The  spatial  eigenvalue  method  is  in  the  class  of  spectral  methods;  it  is 
applied  here  only  to  the  linear  problem  (extension  to  the  nonlinear  case  is 
underway).  Thus,  Kq  must  be  small  in  some  sense.  Although  the  spatial 

eigenvalues  method  was  intended  originally  for  application  to  the  "small"  Re 
range,  say  Re  <  100,  which  is  the  range  emphasized  in  this  report,  in  princi¬ 
ple  it  can  be  applied  to  any  Re;  actually  only  Re  <  1000  has  been  considered. 

The  CPU  time  on  the  VAX  8600  is  of  the  order  of  minutes,  increasing  with  Re; 
for  Re  =  10,  f  =  0.05  and  aspect  ratio  =  3.0,  it  is  10  seconds  if  6  eigen¬ 
values  are  used  (with  a  maximum  of  4  iterations)  and  30  seconds  if  15  eigen¬ 
values  are  used  (with  a  maximum  of  8  iterations).  The  side  moment  coefficients 

computed  with  6  or  15  eigenvalues  differ  by  0.153»  for  this  case.  For  Re  =  1000 

and  f  =  0.1,  the  CPU  time  is  approximately  30  minutes. 

II.  FORMULATION  OF  THE  PROBLEM 

Consider  the  flow  of  a  viscous  fluid  of  kinematic  viscosity  v  and  den¬ 
sity  p  in  a  circular  cylinder  of  radius  a  and  length  2c.  The  aspect  ratio  is 
A  =  c/a.  The  undisturbed  fluid  is  assumed  to  be  in  rigid  body  rotation  with 
angular  velocity  n,  as  observed  in  an  aerobal 1 i Stic  reference  frame.  S';  see 
Reference  6,  e.g.  With  respect  to  the  inertial  reference  frame,  S,  used  here 

the  angular  velocity  is  n  +  t  cos  about  the  x-axis,  where  x  is  the 

dimensional  coning  frequency;  it  is  used  to  define  the  Reynolds  number 

Re  =  (q  +  T  cos  K^)  a2/v  .  (2.1) 

In  S'  the  non-dimensional  coning  frequency  is  i'  =  x/fi  and  the  Reynolds  number 
is  Re'  =  qa2/u.  For  reference,  in  S  the  non-dimensional  coning  frequency  is 
f  =  T  =  t'/(1  +  t'  cos  Kq)  and  Re  =  Re'/(1  -  t  cos  K^).  Using  the  assumption 

already  introduced,  Kq  <<  1,  cos  Kq  is  replaced  by  unity  in  the  above  and  the 

following.  The  typical  time,  length,  velocity,  and  pressure  scales  used  to 

.  - 1 

non-dimensional ize  the  flow  variables  are  (n  +  i)  ,  a,  (n  +  T)a  and 
pa  2  +  x)  2. 
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As  for  any  oscillatory  system  two  problems  are  defined:  the  free  oscil¬ 
lation  or  eigenvalue  problem  and  the  forced  oscillation  problem  here  called 
the  moment  problem.  For  the  eigenvalue  problem,  the  complex  perturbation 
velocities  and  pressure  in  (1.1)  are  written 

{u,v,w,p)  =  U(r)  sin  Kx,  V(r)  sin  Kx,  W(r)  cos  Kx,  P(r)  sin  Kx  .  (2.2) 

From  (1.2)  the  radial  variations  (LI,V,W,P)  satisfy 

[Re’^  (a  -  r“^)  -  iM]  U  +  2  [1  +  (im/r2  Re)]  V  -  P  =  0 
1  ^ 

[Re”^  (a  -  r"^)  -  iM]  V  -  2[1  +  (im/r2  Re)]  U  +  =  0 

1  r 

[Re"^  A  -  iM]  W  -  KP  =  0 
1 


(rU)^  -  imV  -  KrW  =  0  .  (2.3a,b,c,d) 


where  M  =  f  -  m, 


and  A  i  3  +  7  3  -  [—  +  K2]  . 

1  rr  r  r  ^2 

See  also  Reference  8  for  these  equations  in  a  different  context.  The  above 
equations  must  be  solved  subject  to  the  no-slip  condition  atr=l,  U  =  V  =  lfl 
=  0,  together  with  the  appropriate  boundary  conditions  at  r  =  0,  which  depend 
on  m  and  can  be  written 


U  =  V  =  -^  =  0,  m  =  0 


U-iV=W=P=0,  m=l 


U=V=W=0,  m>l.  (2.4a,b,c) 


For  a  discussion  of  the  origin  of  these  boundary  conditions,  see  Batchelor  and 
Gill^  and  Reference  4.  The  system  (2.3),  the  no-slip  condition  at  r  =  1  and 
one  of  (2.4)  constitute  an  eigenvalue  problem  for  K  =  K(m,Re,f).  Since  the 
problem  is  defined  on  a  finite  interval,  it  is  expected  that  the  spectrum  will 
be  infinite  and  discrete.  If  neutral  disturbances  were  possible,  then  K  would, 


'  k  ■■  >1  r 
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of  course,  be  real.  Hie  numerical  solution  of  the  eigenvalue  problem  is  dis¬ 
cussed  in  the  next  section;  here  it  is  assumed  that  an  infinite  sequence  of 
eigenvalues  {K^}  exists  and  can  be  ordered  in  some  sensible  manner. 

For  the  moment  problem,  it  can  be  shown  that  only  m  =  1  is  required.  The 
pivot  point  for  the  coning  motion  is  midway  along  the  axis  of  the  cylinder 
( t  =  0  in  the  notation  of  Reference  4).  The  boundary  conditions  at  r  =  1  and 
X  =  ±A  for  (u,v,w)  were  derived  in  Reference  4;  they  are 


u  =  -i  (1  -  f)  X 
V  =  -(1  -  f)  X 

w  =  i  (1  -  f)  r  .  (2.5a,b,c) 


It  is  convenient  at  this  stage  to  transfer  che  i nhoniogeneous  boundary  condi¬ 
tions  to  the  endwalls  by  subtracting  out  a  suitable  particular  solution  to  the 
perturbed  Nav i er-Stokes  equations.  Let 

u  =  -i  [1  -  f]  X  +  u  (r,x) 

y  =  -[1  -  f]  X  +  V  (r,x) 

w  =  i  [1  -  f]  r  +  o(r)  +  w  {r,x) 

g  =  -[1  -  f2j  rx  +  p  (r,x) 

J 1  ( Ar ) 

a(r)  =  2if  [r  -  — - ]  (2 . 6a  ,b ,c  ,d  ,e) 

J^(  a) 

where  is  the  Bessel  function  of  the  first  kind  and  order  1  and 

A  =  (1  +  i)  [(1  -  f)  Re/2]^^^.  The  boundary  conditions  at  r  =  0,1  and  x  =  ±A 
in  terms  of  (u,v,w,p)  become 


u-iv=w=p=0,  r=0 

u=v=w=0,  r=l 


u 


V 


0,  w  =  -o(r) 


X  =  i-A 


(2. 7a  ,b,c ) 


and  the  differential  equations  satisfied  by  (u,v,w,p)  are  (1.2)  with  m  =  1: 


i  [f  -  1]  u  -  2v  =  -p  +  —  |72u  -  I  , 

r  Re  1  r-2  ,.2  f 


r2 


i  [f  .  1]  V  +  2u  =  +  1-  /yZv  .  .  llu  i 

r  Re  '  r2  r^  ‘ 


i  [f  -  1]  w  =  -p  +  —  /  72  w  -  15L  I 
^  Re  '  r2) 


r2 


(ru)^  -  iv  +  rw^  =  0 


{2.8a,b,c,d) 


Assuming  that  (u,v,p)  =  (U(r),  V{r),  P(r))  sin  Kx  and  w  =  W(r)  cos  Kx,  it 
follows  that  U,V,W,P  satisfy  (2.3)  with  m  =  1;  the  odd  and  even  variations 
with  X  follow  from  (2.5)  and  (2.6).  Let  the  infinite  sequence  of  eigenvalues, 
{K^},  with  m  =  1  be  denoted  by  {k^}  and  the  eigenfunctions  by  u^.  v^,,  w^, 

Pf,.  The  solution  of  (2.8)  is  expressed  as  an  eigenfunction  expansion: 


n 


n 


a  sin  k  X 

X)  - p  (r)  .  (2.9a,b,c,d) 

=  1  sin  k  A 
n 


w  =  - w  (r)  ,  p  = 

n  =  1  sin  k  A 
n 


o  sin  k  X 
n _ n_ 

1  sin  k  A 
n 


Vn(r), 


u  =  X 
n  =  1 


a  s  i  n  k  X 
_n _ ri_ 

sink  A 
n 


u„{r), 


V  = 


Discussions  of  how  {k^ }  can  be  sensibly  ordered  and  how  the  as  yet  undeter¬ 
mined  constants  {a^ }  can  be  found,  by  applying  the  endwall  boundary  condi¬ 
tions,  are  given  in  the  next  section. 


III.  CALCULATION  OF  THE  EIGENVALUES  AND  SERIES  COEFFICIENTS 

The  first  step  in  the  calculations  is  to  find  {k  }  and  {u  ,v  ,w  ,p  },  the 

^n  ^nnn^n 

eigenvalues  and  eigenfunctions  of  (2.3)  with  m  =  1,  subject  to 


TV  =  W  =  p 

n  n  ^n 


=  0, 


1  =  V  =  w  =  0, 


r  =  0, 

r  =  1. 


(3. la  ,b ) 


Two  independent  methods  were  used  to  generate  the  eigenvalues;  since  both 
iterative  methods  are  well  documented  elsewhere,  only  a  brief  description  of 
them  is  given. 

Firstly,  the  complete  orthonormalization  procedure  of  Davey,iO  as 
implemented  by  Kitchens,  Gerber,  and  Sedneyii  was  employed.  The  numerical 
integration  was  carried  out  using  a  fourth  order  Runge-Kutta  scheme.  Secondly, 
the  compact  finite  difference  scheme  of  Malik,  Chuang,  and  Houssaini^^  was 
used.  In  both  methods  the  singularity  in  the  equations  at  r  =  0  was  removed 

by  transferring  the  boundary  conditions  to  r  =  e  with  0  <  e  <<  1  and  using  the 

Taylor  series  solution  given  in  Reference  11.  The  methods  gave  consistent 
results.  The  complete  orthonormalization  method  was  more  efficient  at 
Reynolds  number  of  order  102.  /\t  higher  Reynolds  nunbers  where  some  of  the 

eigenfunctions  begin  to  vary  rapidly,  the  finite  difference  method  was  found 
to  be  the  more  efficient  method. 

With  both  methods  a  sufficiently  accurate  guess  for  kp  was  required  if 

the  iterations  were  to  converge.  This  was  a  hindrance  initially  since  little 

previous  knowledge  of  the  distribution  of  these  eigenvalues  was  available;  in 
certain  cases  it  was  necessary  to  calculate  a  large  number  of  the  sequence 
{k^}.  To  obtain  first  guesses  techniques  were  developed  to  determine  approxi¬ 
mations  to  the  eigenvalue  distribution  in  certain  special  limits,  and  if 
necessary,  follow  their  evolution  into  the  required  regime. 


1.  A  RELATED  TAYLOR  VORTEX  PROBLEM. 

Firstly,  previous  knowledge  of  the  eigenvalue  distribution  for  the 

classical  Taylor  vortex  problem  was  used  (see  Reference  5).  Consider  a  basic 
flow  driven  by  the  rotation  of  the  cylinder  r  =  e  while  the  outer  cylinder  is 
held  fixed.  No-slip  boundary  conditions  are  applied  at  r  =  e  and  at  r  =  1. 
The  frequency,  f,  and  azimuthal  wave  numbers,  m,  are  now  set  to  zero;  at 

sufficiently  small  supercritical  Reynolds  numbers  the  resulting  eigenvalue 

distribution  is  as  sketched  in  Figure  1.  The  two  real  wave  numbers  k^  and  k2 

correspond  to  the  two  neutrally  stable  perturbations  which  exist  in  this 

regime.  In  addition,  there  are  the  following  three  sets  of  eigenvalues: 


(1)  A  pure  imaginary  set  {kj^^}. 


n  =  1,2,.. 


(2)  A  complex  set  of  eigenvalues 
with  k3p+2  “  ^  3n+l’ 


(3)  A  complex  set  of  eigenvalues  k3p+2  , 


n  —  1,2,... 


being  in  the  first  quadrant  with  real  part  -  n. 

The  eigenvalues  and  k2  can  be  associated  with  (3)  and  (2),  respectively, 

and  at  subcritical  Reynolds  numbers  move  into  the  complex  plane  with  k^  =  C  2* 


At  higher  Reynolds  numbers  pairs  of  eigenvalues  (^3^+1,  k3n+2)  '^erge  on  the 

real  axis  and  then  move  along  the  real  axis.  Blennerhasset  and  Hall  5  found 
that  the  eigenvalues  {^30+1,  k3p+2,  ^3^+3  }  for  n  =  0,1,2,  have  azimuthal 

velocity  eigenfunctions  with  n  zeros  in  (e,l).  Thus,  in  the  axisymmetric 
problem,  there  is  a  rational  method  of  ordering  the  eigenvalues,  in  triplets. 
If  m  or  f  is  now  taken  to  be  nonzero,  the  symmetry  of  the  eigenvalue  distri¬ 
bution  described  above  is  broken  but  the  ordering  from  the  axisymmetric  steady 
Taylor  vortex  problem  can  be  retained  while: 

a.  f  and  m  are  gradually  increased  from  zero  to  their  required  values. 

b.  The  basic  flow  is  gradually  varied  from  that  appropriate  to  the 
Taylor  problem  to  rigid  body  rotation. 

c.  The  boundary  conditions  at  r  =  e  are  gradually  varied  from  no-slip  to 
those  given  by  the  Taylor  series  expansion  of  Reference  11. 

It  was  found  that  it  was  not  necessary  to  perform  a,  b,  and  c  independ¬ 
ently  and  that  having  found  a  set  of  eigenvalues  {k  }  for  some  f  and  Re,  other 

sets  could  be  generated  by  varying  f  and/or  Re  slowly,  i.e.,  by  extrapolation. 
However,  at  values  of  Re  greater  than  about  100,  the  eigenvalues  (2)  and  (3) 
were  found  to  be  very  sensitive  to  changes  in  Re  and  it  proved  rather  ineffi¬ 
cient  to  calculate  them  by  this  method.  For  large  Re  an  asymptotic  approach 
was  used,  described  below.  Another  asymptotic  approach,  for  small  Re,  is 
necessary  and  outlined  next. 

2.  Re  +  0. 

For  values  of  Re  <  10,  approximately,  estimates  for  the  eigenvalues  can  be 

found  by  solving  the  Re  +  0  problem,  i.e.,  the  Stokes  limit.  Essentially,  the 
same  procedure  is  used  as  for  the  classical  Stokes  flow  solution  but  Re  0  is 
a  regular  perturbation  for  the  internal  flow  under  consideration.  In  addition 
to  the  momentum  and  continuity  equations,  (2.8),  a  derived  equation  is  used. 
The  linearized  Na vi er-Stokes  momentum  equations,  in  vector  form  with  velocity 
vector  q*  =  {u*,v*,w*),  see  (1.1),  can  be  written 


i (f  -  1)  q*  +  2  (i^  X  q* )  =  -grad  p*  -  (1/Re)  curl 

where  oo*  is  vorticity  and  i^  is  the  unit  vector  in  the  x  direction.  Using 

div  q*  =  div  curl  u*  =  0  and  the  definitions  (1.1)  yields 


r  [v2  p  -  p/r2]  =  2  [iu  +  (rv)^]  .  (3.2) 

Next,  the  scaled  pressure  P  =  g  Re  is  introduced  and  it  is  convenient  to 
change  variables  u  and  v  to  U  =  u  +  iv  and  V  =  u  -  iv.  A  regular  perturbation 


(U.V.w.P)  =  E  (U..  V..  w..  P.)  ReJ 
j  =  0  J  J  J  J 


is  substituted  into  (2.8),  (3.2),  and  (2.7)  and  the  latter  is  expanded  in  a 
power  series  in  Re.  The  result  is 


,2  U  ,  .  +  Pj/r 


,2  Vj  -  (4/r2)  Vj  .  Pj,  -  Pj/r 


72  .  (l/r2)  .  P,^ 


Uj,  *  »jr  *  (2/'-)  *  2  "Jx  '  “ 


v2  P  -  (l/r2)  P  =  0  . 

J  J 


(3.3a,b,c,d,e) 


Using  the  expanded  boundary  conditions  (2.7),  the  trivial  solution  is 
obtained  for  the  0(1)  terms.  The  0(Re)  solution  is  determined  from 

homogeneous  and  non-homogeneous  solutions  of  (3.3),  with  respect  to  P, . 
Expressing  the  0(Re)  terms  as  a  modal  decomposition,  as  in  (2.2), 

(Uj ,V^ ,w^ ,P^)  =  F^(r)  sin  x^x,  G^(r)  sin  x^x,  H^(r)  cos  X^x,  n^(r)  sin  x^x  , 

the  eigenfunctions  are  found  to  be 


^n  ^  '^1  ^^n^  '^0 

Op  =  [f  ^2  ^^n^  ^^n^  '^2  ^^n^ 

H„  =  f  [r  Jj  (A„)  J(,  (»„r)  -  Jj  (A„)  Jj  (A„r)]/Jj  (A„) 

"n  =  2  Oj  (A„r) 

where  =  ix^,  are  the  roots  of 

Ap  Jf  (Jq  -  ^2^  ■  ^  ^n  *^2  ■  ^ 


(3.4a,b,c ,d) 


(3.5) 


V“’ 

v*' 

A  v.vNr"'* 


Obviously,  the  x  do  not  depend  on  any  of  the  parameters  of  the  problem;  the 
n 

first  several  x  are  shown  in  Figure  2,  indicated  by  k-  since  X  =  lim  k  . 

n  ^  •'n  n„.  „n 


These  eigenvalues  separate  into  three  branches  or  sets: 

(1)  a  set  which  is  pure  imaginary, 

(2)  a  set  in  the  fourth  quadrant,  and 

(3)  a  set  which  is  conjugate  to  (2). 


Re+0 


As  an  alternate  notation  for  the  x's  or  k's,  let  the  eigenvalues  be 


denoted  by  x  where  £=  1,2,3  is  the  branch  and  n  =  1,2,...  is  the  index 
t,n 


along  a  branch.  If  Jq  ^  and  ^2  n  zeroes  of  Jq  and  J2, 

respectively. 


'l.n 


^  ^'^0,n+2  '^2,n^^^ 


for  n  >>  1.  This  asymptotic  estimate  is  quite  accurate  for  n  >  1.  From 

this  x^  ■  ^1  n  “  ’  ^tually,  the  eigenvalues  for  finite  Re  satisfy 

Im  (k  ,  -  k  )  •>  ±  ir,  +  for  i  =  1,3  and  -  for  i  =  2,  for  all  i  and  Re  if 
£,n  +  l  i,n' 

n  >  nQ(Re).  The  relationship  nQ{Re)  is  only  known  empirically.  An  algorithm 
for  obtaining  a  first  guess  for  k^  ^  is: 

For  given  i.  Real  k  is  obtained  by  extrapolating  Real  k  _  and 

Real  and  Im  ±  ,  for  n  >  n^. 

The  k  were  computed  up  to  Re  =  100,  starting  from  the  X  ,  by  slowly 
t ,n  £  ,n 

increasing  Re  and  using  the  algorithm.  Relatively  large  increments  in  Re 

could  be  used  for  i  =  1,  but  the  process  became  tedious  for  £  =  2,3.  Results 

for  Re  =  100,  f  =  0.1,  and  n  <  5  are  shown  in  Figure  2.  The  first  eigenvalue 

on  each  branch  and  k3^2  labeled.  Note  the  separation  of  k3^2  *^3,1* 

increases  monotonically  with  Re;  lk3^2  "  *^3,ll/l^3  2  "  ^3  1  ^  ” 

10  and  4.62  for  Re  =  100,  both  for  f  =  0.1.  The  ordering  of  the  k  is 

A 

determined  by  the  numoer  of  zeroes  of  the  eigenfunctions:  n  zeroes  corre¬ 
sponding  to  ^  1  zeroes  corresponding  to  k2^n  *^3,n’ 

ordering  is  consistent  witn  the  analytical  results  obtained  in  the  Stokes 


In  calculating  the  flow  variables,  the  must  be  taken  in  groups  of 

three:  i  =  1,2,3  for  n  =  1,  n  =  2,  etc.  Departure  from  this  ordering 

introduces  errors;  the  error  depends  mainly  on  the  degree  of  departure,  n  and 

on  the  method  used  to  calculate  a  . 

n 

A  more  efficient  method  of  estimating  k2^n  *^3,n  large  Re  was 
needed  and  for  this  purpose  an  asymptotic  approximation  was  derived. 


An  asynptotic  solution  to  (2.3),  with  m  =  1  and  K  =  k,  can  be  found  by 
expanding  the  disturbance  field  and  wave  numbers  in  the  form 


k  =  >/te  ICq  +  +  ...  , 


^  "  “o  Re  ’ 


=  ''o  ^  “•  • 


1.1 

w  =  - w.  + - w,  +  ...  , 

U  O  /  O  A 


Re 


3/2 


=  l-p  +J_ 

Re  ^  Re  2 


P,  +  ••• 


(3.6a,b,c,d,e) 


The  zeroth  order  approximation  to  (2.3a,b)  then  yields 

{-<1  -  if  +  i}  Uq  +  2vq  =  0 

-  if  +  i  }  Vq  -  2uq  =  0  (3.7a,b) 


and  the  consistency  of  (3.7a,b)  requires  that  either 


•**a 


•&I 


a 


where  Real  Kq  >  0  is  required  so  that  (3.8a)  and  (3.8b)  correspond  to  branches 
2  and  3,  respectively.  A  spectrum  is  not  obtained  in  this  approximation  to  k. 

Corresponding  to  (3.8a),  the  zeroth  approximation  to  (2.3c,d)  yields 


*^0  "" 


(3.9a,b) 


where  '  =  d/dr;  uq  is  obtained  as  follows. 


The  next  order  approximation  to  (2.3a,b)  yields  a  pair  of  inhomogeneous 
linear  equations  for  u^,  v^.  The  operators  in  the  homogeneous  parts  of  these 

equations  are  identical  to  those  in  the  l.h.s.  of  (3.7);  thus,  the  inhomogene¬ 
ous  terms  must  satisfy  a  solvability  condition,  which,  for  k,.  from  (3.8a), 
reduces  to  ^ 


u^  +  (1/r)  u^  +  Uq  =  0 


(3.10) 


where  =  -2<g<^  (l+f)/(2+f)  and  the  required  solution  is 


Uq  =  Jg(Dr). 


The  <Q  from  (3.8b)  leads  to 


Uq  +  (1/r)  Uq  +  (B^r^  -  4)  Ug/r^  =  0, 


(3.11) 


where  =  '2<gic^  (3-f)/(4-f)  and  the  required  solution  is 


Uq  =  02(Br), 


(3.12) 


The  boundary  condition  uq(1)  =  0  determines  for  each  branch;  their  values, 
corresponding  to  (3.8a,b,)  are  given  by 

=  -  (v^/4)  (Hi)(2+f)  /  (l  +  f)^/2 


=  .(^/T/4)(l-i)(4.f)  /  (3-f)' 


( 3. 13a  ,b) 


V 


>1 


w.vL'vw- n* 


j»  • 


» 


respectively,  s  =  1,2,...  .  Thus  (3.6a),  using  (3.8a)  and  (3.13a)  or  p.8b) 
and  (3.13b),  approximates  the  eigenvalues  on  branches  2  and  3,  respectively, 

correct  to  0(Re"^/^).  From  (3.9b)  and  the  boundary  condition  Wq(1)  =  0,  a 

boundary  layer  of  thickness  0(Re"'^/^)  must  be  introduced  at  r  =  1  but  this 

does  not  affect  given  by  (3.13).  Further  asymptotic  anproximations  to  the 

eigenvalues  will  be  reported  in  a  separate  publication. 

Consider  a  numerical  example:  Re  =  1,000,  f  =  0.1  and  branch  3.  The 

first  term  of  (3.6a),  with  (3.8b),  gives  V  1,000  kq  =  38.0789  (1+i)  while 

the  first  two  terms,  using  (3.8b)  and  (3.13b)  with  s  =  1,  give  an  approxi¬ 
mation  37.8460  +  i  38.3117.  Using  these  as  first  guesses  in  the  eigenvalue 
iteration  routine,  the  iteration  process  diverged  for  the  former  but  converged 
for  the  latter,  even  if  rounded  to  2  decimal  places.  The  converged  value  is 
37.8471  +  i  38.3196.  The  iteration  process  also  converged  using  guesses 
obtained  from  the  first  two  terms  of  (3.6a)  and  (3.13)  with  s  =  2  and  3  but 
not  with  s  =  4.  If  the  third  term  in  (3.6a)  is  also  included,  but  the 
boundary  layer  at  r  =  1  neglected,  the  iteration  process  converges  for  s  <  5. 

Although  the  first  term  in  (3.6a)  does  not  provide  a  first  guess  which 
converges  to  an  eigenvalue,  it  can  be  used  to  obtain  a  "differential  correc¬ 
tion"  to  an  eigenvalue.  In  fact,  an  extrapolation  of  the  eigenvalue  for  Rej, 

f]'  to  Re2,  f^  can  be  derived  from  the  first  term  in  (3.6a)  which  is  second 
order  correct  in  ARe  =  Re2  -  Re^  and  Af  =  f2  -  f],.  For  branch  2,  the 
increment  in  k  is 


Ak 


=  (l-i )/2>/r  [aRo  +  Re^  Af/d+fp)] 


and  for  branch  3 


Ak 


=  [V3^  (l+i)/2  VT  >/^o]  [aRo  -  Rep  Af/(3-fo)] 
where  Reg  and  fg  are  the  average  values  of  Re2,  Re;^,  and  f2,  f^,  respectively. 


Having  calculated  just  three  eigenvalues  along  each  of  the  branches  2  and 
3,  the  first  guess  for  the  next  one  was  obtained  using  the  asymptotic  approxi¬ 
mation  and  extrapolation  in  n.  For  larger  n,  the  algorithm  stated  previously 
was  used.  Although  the  asymptotic  method  for  Re  ♦  •  described  here  does  not 
give  approximations  for  branch  1,  those  eigenvalues  can  be  obtained  rather 
easily  by  extrapolation  from  Re  =  100.  In  Figure  3,  the  converged 
eigenvalues,  k  ,  are  presented  for  Re  =  1,000,  f  =  0.1,  i  -  1,2,3  and 
n  <  24. 


4.  CALCULATION  OF  o^. 

Reverting  to  the  first  notation  for  the  eigenvalues,  assume  that  {k^}  are 

known  for  1  <  n  <  3M,  where  M  is  a  positive  integer.  Truncating  the  series  in 

(2.9)  at  n  =  3M,  the  coefficients  {o^}»  1  <  n  <  3M,  were  determined  by  either 

of  two  ways.  The  first  is  collocation.  Let  {r.},  1  <  j  <  M,  be  a  sequence  of 
points  such  that  ^ 

0  <  <  r2  . . .  <  <  1. 

Enforcing  the  boundary  conditions  (2.7)  at  the  points  rj^,  ....  r|«|  leads  to 

3M 

‘'■j>  =  Z  “n  ('•jl  -  “ 

n  =  1 
3M 

^2  <''j>  ^  Z  “n  'n  <''j>  ■  “ 
n  =  1 

3M 

E3  (Cj)  =  ^  (a^  cot  k^A)  w^  (r^)  +  o(rj)  =  0  .  (3.14a,b,c) 

n  =  1 

These  31‘1  linear  equations  are  solved  for  the  coefficients  ap  ...  The 

success  of  the  method  requires  that  the  flow  properties  converge  in  some  sense 
as  M  increases.  Experience  has  shown,  for  Re  <  1,000,  that  the  required  value 

of  M,  to  achieve  a  yiven  accuracy,  increased  monotoni cal ly  with  Re.  Numerical 
examples  are  given  in  Section  V. 

The  second  is  a  least  squares  method.  The  "normal"  equations  for  Op  are 
obtained  by  minimizing  the  error 

g(an)  =  {|Ei(r)l2  +  |E2(r)l^  +  |E3(r)l2|r®  dr  .  (3.15) 

^  0 

Exponents  e  =  0,1  gave  essentially  the  same  results.  A  measure  of  the  error, 
E,  in  Op  is  given  by 

1 

=  g/  (  1 o( r ) I ^  r^  dr.  (3.16) 

•^0 

Tne  coefficients  of  the  linear  equations  for  are  given  in  the  Appendix. 


IV.  PRESSURE  AND  MOMENT  COEFFICIENTS 


The  velocity,  pressure  and  moment  exerted  by  the  liquid  on  the  container 
are  quantities  of  physical  interest  calculated  by  the  spatial  eigenvalue 
method.  The  last  is  of  direct  use  in  the  application  of  the  theory  to  the 
study  of  projectile  stability;  e.g.,  it  can  be  used  to  calculate  yaw  growth 
rate. ‘♦’13  Measurements  of  pressure  have  been  made  mostly  on  the  endwall  and 
the  sidewall,  an  exception  being  those  on  the  axis  of  the  cylinder  by 
Aldridge,  i"*  Only  the  endwall  case  is  considered.  The  pressure  transducer 
measurements  are  processedis  to  produce  amplitude  and  phase.  This  amplitude 
is  compared  with  calculations  in  the  next  section.  Comparisons  are  made  using 
only  non-dimensional  pressure  and  moment  coefficients.  In  Reference  16 
dimensional  quantities  are  used;  these  are  converted  to  dimensionless  forms 
for  present  purposes.  Derivations  of  pressure  and  moment  coefficients  are 
outlined;  details  can  be  found  in  References  4  and  13. 

The  pressure  calculated  from  the  Navier- Stokes  equations,  pj^^,  is 

Pf^3  =  1/2  r^  -  p*  (4.1) 

in  the  linear  theory.  However,  since  the  measurement  is  made  in  a  reference 

frame  fixed  to  the  cylinder,  with  coordinates  r,e,x  ,  the  coordinates  used  in 
the  analysis,  r,0,x,  must  be  transformed  to  this  frame.  For  example, 

1/2  r^  =  1/2  r^  -  Kq  r  X  cos  (ft-e)  +  0  (K^)  .  (4.2) 

The  disturbance  pressure,  Ap ,  is  the  pressure  not  including  the  contribution 
from  solid  body  rotation.  At  points  fixed  on  the  surface  of  the  cylinder 

Ap  =  1/2  r^  -  p*  (r,e,x,t)  -  1/2  r^  . 

Using  p*  (r,d,x,t)  =  p*  (r,0,x,t)  +  0(Kq),  (1.1)  and  (4.2), 


Ap  =  -K^  [-pj  sin  (ft-e)  +  (P|^  +  rx)  cos  (ft-e)]  (4.3) 


where  sub  R  and  I  denote  real  and  imaginary  parts  of  p,  respectively,  and 

the  rx  factor  in  (4.2)  can  be  replaced  by  rx  without  changing  the  order  of 
approximation;  (2.6d)  relates  pj  to  P[^,  pj  and  the  latter  are  obtained 

from  (2.9d).  Note  that  p^^  +  rx  =  P[^  +  f2rx  and  Pj  =  pj.  The  pressure 

coefficient,  Cp,  is  defined  as  the  amplitude  of  Ap/K^,  i.e., 

r  r/  .  2.1/2 

Cp  =  [(Er  +  rx)  +  Pj  ]  . 
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(4.4) 


Certain  limiting  forms  of  Cp  can  be  derived.  Firstly,  consider  f  ♦  0. 
From  (2.7c),  the  boundary  condition  is 


w  =  fu  (r)  ,  X  =  ±A 


so  that  the  from  (3.14)  and  (u,v,w,p)  in  (2.9)  are  0(f).  Therefore, 

(u,v,w,p)  =  f  (u°,v°,w°,p°)  +  O(f^) 


where  (u°,  v° ,  w° ,  p°)  expand  as  in  (2.9).  Tlie  eigenvalues  are  evaluated  for 
f  =  0  and  the  endwall  boundary  condition  is 


=  u  (r)  ,  X  =  ±A  , 


which  leads  to 


Cp  ■  lf| 


(4.5) 


For  Re  =  8.8,  A  =  3.148  and  r  =  0.667,  (4.5)  gives 


Cp  =  .82  Ifl  +  O(f^) 


For  large  values  of  f,  the  endwall  boundary  condition  becomes 


w  =  -2ifr  +  0  )  , 


X  =  ±A 


for  r  >>  Ifr^^^;  this  is  obtained  by  using  the  asynptotic  expansions  for 
Ji ( xr)  and  Ji ( x)  in  (2.7c).  Thus,  the  flow  has  a  boundary  layer  structure  of 

thickness  0  (ifj  )  at  r  =  0  which  implies  that  e  (see  Section  III)  must 

be  chosen  such  that  e  <<  |f("^/^.  From  (3.14c)  and  (2.9d)  p  ~  |fj  and  a  first 

approximation,  neglecting  the  0  dfl"^^^)  boundary  layer,  gives 


Cp  =  f"  r  A  +  0  (Ifl) 


(4.6) 


on  the  endwal  1  . 


•  ^  *  *  m  •  »  •  M  ^  ?  «.  »  -V  «  •  -  •  ^  ^  •  fcN  _  t  fcN  «  »  .Ni 
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If  f  =  1,  the  perturbation  is  zero  which  is  clear  from  (2.6),  (2.7),  and 
(2.8).  Therefore,  (4.4)  gives 


Cp  =  rx  ,  f  =  1  , 


(4.7) 


a  linear  variation  on  the  endwall  and  sidewall.  Since  (4.7)  is  exact,  it  can 
be  used  in  error  analyses  of  the  numerical  results. 


For  r  <<  1,  it  can  be  shown  that,  on  the  endwall,  for  all  Re  and  f 


Cp  =  0(r). 


(4.8) 


The  moment  exerted  by  the  liquid  on  the  container  can  be  calculated  from 
the  rate  of  change  of  angular  momentum^* ^  or  from  the  integration  over  the 
cylinder  surface  of  the  pressure  and  shear  forces. *♦’ 16*  17  jhg  latter  is 
outlined.  The  three  components  of  the  moment  are  evaluated  in  an  aeroballis- 
tic  coordinate  system.  The  axial  or  roll  component  is  zero  to  0(Kq);  the 

transverse  components  are  0  (Kq)  and,  to  that  order,  are  the  same  in  the 

aerobal 1 istic  and  inertial  frames.  The  transverse  components  can  be  separated 
into  the  overturning  or  side  moment  and  the  in-plane  moment  which  acts  to 
change  the  coning  rate.i^  The  side  moment  is  important  for  projectile 
stability,  as  shown  by  Murphy, and  is  the  only  component  considered.  The 
following  definitions  are  used;  nondimens ional  side  moment  coefficient  = 

^LSM  =  moment/2iTp  a  c  (t  +  Ji)  ^  Ci_5f^  is  expressed  as  the  sum  of 

the  four  contributions;  the  pressure  on  the  sidewall,  press  ire  on  the 
endwalls,  shear  force  on  the  sidewall  and  shear  force  on  the  endwalls.  These 
are  given  by 


V  X  gj  (1  ,x)  dx 

J  n 


-(fA)' 


i 

'  5  El 


A)  r^  dr 


•(fARe) 


A 

■u 


3  (xvj^  -  Wj)/3r 


+  1  -  f 
r=l 


fRe)  Real  l-f  +  \  3(v-iu)/3x 


(4.9a ,b ,c ,d ) 


respectively.*  For  the  reasons  given  in  deriving  (4.7)  C|_5f/|  =  0  for  f  =  1. 

The  derivative  in  (4.9d)  will  introduce  a  factor  kp  in  the  series  used  to 

calculate  (4.9d)  which,  it  was  conjectured,  would  adversely  affect  conver¬ 
gence.  By  integrating  (2.8a,b)  over  x  and  r,  an  alternative  form  was  derived 
which  did  not  have  this  factor.  It  was  found,  numerically,  that  the  two  forms 
agree  to  2%  or  less  for  |kp|  ~  100. 

If  the  pressure  and  moment  coefficients  are  defined  using  the  spin,  n, 
and  non-dimensional  coning  frequency,  t',  relative  to  S',  and  denoted  by  Cp' 

and  Clsm'.  they  are  related  to  Cp  and  Clsm  by  Cp  =  Cp*  (1  -  f  cos  Ko)2  and 

^LSM  =  ^LSM'  (1  -  ■f'  (^os  Kq). 

V.  RESULTS 

A  sampling  of  results  obtained  by  the  spatial  eigenvalue  method  is 
given.  Some  of  these  are  used  to  discuss  the  numerical  process  but  most  of 
the  results  presented  either  have  intrinsic  interest  or  are  used  to  compare 
with  finite  difference  calculations  or  experiment.  A  standard  set  of  the 
basic  parameters  was  chosen  to  facilitate  comparison  of  several  of  the 
figures;  these  are  Re  =  10,  f  =  0.1,  A  =  3.0.  Since  the  original  motivation 
for  this  work  was  the  calculation  of  Cp  and  C|_5[.i  for  "small"  Re,  this  range  is 

emphasized  here.  For  "large"  Re,  the  theories  of  References  4,  and  13  are 
aval labl e. 

Instead  of  an  analysis  of  the  numerical  processes  used  in  the  spatial 
eigenvalue  method,  a  number  of  empirical  studies  were  made,  and  some  con¬ 
clusions  based  on  these  are  presented.  Recall  that  two  methods  were  used  for 
determining  the  eigenvalues,  viz.,  orthonormal ization i O’ i i  and  finite 
di fference, 12  and  that  two  methods  were  used  to  determine  the  coefficients, 
a^,  viz.,  collocation  and  least  squares.  As  expected,  results  were  not  very 

sensitive  to  the  choice  of  method  for  computing  k  .  Extensive  calculations 

t,n 

were  made  with  a  program  combining  orthonormalization  with  least  squares  and 
another  combining  the  finite  difference  technique  with  collocation;  these  are 
denoted  by  LS  and  COL,  respectively. 

Of  the  physical  parameters,  the  calculation  is  more  sensitive  to  Re  and  f 
than  A;  recall  that  the  k^  ^  are  independent  of  A.  Of  the  numerical  parame¬ 
ters,  obviously  the  number  of  eigenvalues,  N  =  3M ,  is  most  important;  in  COL 
the  number  of  collocation  points  is  M.  In  all  calculations  using  collocation, 
the  collocation  points  were  equally  spaced  (exactly  or  to  within  ±1  of  the 
point  number).  The  number  of  points  required  in  the  numerical  integration  of 
the  differential  equations  for  the  eigenfunctions  is  directly  related  to  N.  A 
number  of  points  equal  to  7M  was  found  to  be  adequate  although  a  larger  number 
was  usually  employed;  recall  that  the  largest  number  of  zeroes  of  the 
eigenfunctions  is  M.  There  are  several  other  numerical  parameters  having  to 
do  with  the  number  of  terms  in  series,  convergence  of  iteration  processes, 
etc,  which  will  not  be  discussed. 

*Ir.  (4.?,  the  aoKelete  expressions  for  shear  stress  are  used  whereas  approxi- 
natioKS  for  large  Re  are  introduced  in  Reference  13. 


The  first  illustration  is  for  the  calculation  of  C|_5|^  by  LS  and  COL  for 

the  case  Re  =  21.5,  f  =  .0621,  A  =  1.042  and  3  <  N  <  42,  shown  in  Figure  4; 

there  is  a  break  in  the  ordinate  scale  to  accommodate  the  COL  point  for  N  =  9. 
The  LS  results  are  less  sensitive  to  decreasing  N  than  the  COL  results.  The 
LS  value  for  N  =  6  is  within  1.7X  of  the  result  for  large  N  and  required  only 
6  seconds  of  CPU  time  on  the  VAX  8600.  The  LS  and  COL  results  are  essentially 
the  same  for  N  >  12.  A  comparison  with  an  experimental  result  is  shown  later 

for  the  same  Re,  f,  and  A  used  here. 

The  calculated  values  of  Cp  (.667)  vs  N,  using  LS  and  COL,  are  shown  in 

Figures  5  and  6  for  Re  =  100  and  1,000,  respectively,  and  f  =  0.1,  A  =  3.  In 
Figure  5  the  differences  in  Cp  for  N  >  36  are  negligible.  At  N  =  12,  the  LS 

value  is  0.13X  below  and  the  COL  value  is  Al%  below  the  result  for  N  =  72. 
In  Figure  6  the  approach  to  a  limiting  value  for  large  N  is  not  as  clear.  The 
LS  and  COL  results  differ  by  1%  at  most  for  N  >  48.  The  limiting  value  is 

approximately  Cp  (.667)  =  .205.  At  N  =  12  the  COL  value  is  2.4^  below  the 

limiting  value,  and  the  LS  value,  offscale,  is  7.8%  below  the  limiting  value. 

In  all  the  cases  discussed  hereafter,  N  was  chosen  large  enough  to  give 
limiting  values  of  Cp  and  Clsm  sense  used  above,  and  the  values 

obtained  using  LS  and  COL  agreed  to  at  least  three  significant  figures  with 
the  exception  of  the  Cl5m  A  results  for  Re  =  1000.  When  quoting  results 

from  finite  difference  methods,  for  comparisons,  Kq  must  be  specified;  results 

from  the  linear  theory  used  here  are  independent  of  Kq,  of  course. 

The  variation  of  Cp  with  r,  on  the  endwall,  is  shown  in  Figure  7  for  Re  = 

1.0,  f  =  0.1,  and  A  =  3.0.  In  addition,  results'^  from  Strikwerda's  method  are 
presented  for  =  2°.  The  former  is  essentially  linear  with  r  for  approxi¬ 
mately  0  <  r  <  0.2  but  the  latter  is  not.  From  (4.8),  Cp  =  0(r).  This 

comparison  and  some  others  indicate  a  loss  in  accuracy  of  Strikwerda's  method 
for  Re  <  1. 

In  Figure  8,  the  spatial  eigenvalue  calculation  of  Cp  on  the  endwall  vs  f 

is  presented  for  Re  =  10  and  A  =  3.0  at  r  =  0.5  and  1.0.  Recall  that 

Cp  =  0(|f|)  for  f  0,  from  (4.5),  but  the  dependence  is  not  linear  over 

0.1  <  f  <  1.0.  Values  of  Cp  from  the  asymptotic  approximation  (4.6)  are  also 

given.  At  f  =  0.8  and  r  =  1,  the  estimate  (4.6)  is  8%  lower  than  the  calcu¬ 
lation;  at  f  =  1.1  (not  shown)  and  r  =  1,  it  is  3%  higher. 

The  variation  of  Cp  on  the  endwall,  x  =  A,  at  r  =  .667,  over  the  range 

1  <  Re  <  1000  and  for  f  =  0.1,  A  =  3  is  shown  in  Figure  9.  Over  that  range, 

Cp  varies  by  a  factor  of  2.3. 
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The  variation  of  Cp  at  x  =  A,  r  =  .667,  over  the  range  0.5  <  A  <  15  and 

for  Re  =  10,  f  =  0.1  is  shown  in  Figure  10.  In  projectile  application  usually 

I  <  A  <  5;  over  that  range  Cp  increases  by  a  factor  of  1.3. 

A  comparison  of  calculated  pressure  coefficient  vs  f  with  some  experi¬ 
mental  measurements  of  Hepner,  Kendall,  Davis,  and  Tenly^^  is  shown  in  Figure 

II  for  Re  =  3.1  and  A  =  3.148.  The  pressure  transducer  was  located  at  r  = 
0.667  and  the  coning  angle  was  Kq  =  2°.  The  radius  of  the  cylinder  was  3.18 

cm;  the  fluid  had  a  density  of  0.969  gm/cm^  and  a  kinematic  viscosity  of 
60,000  cs.  According  to  Nusca  and  D'Amico^®  the  measured  Cp  and  f  are  in  the 

inertial  reference  frame,  S.  The  data  include  estimates  of  errors  in  Cp,  as 

presented  by  Hepner  et  al.^^  pq^.  almost  all  points,  the  calculated  Cp  is  less 

than  the  measured  Cp;  for  the  larger  f  it  is  outside  the  error  bar.  The 

reason  for  the  bias  in  this  comparison  is  not'  known. 

The  remainder  of  this  section  is  devoted  to  presentation  of  results  of 
calculating  Cls^^  and  some  comparisons.  There  are  more  publications  that 

report  calculation  of  C|_3|^  than  Cp,  starting  with  the  first  application  of  a 

finite  difference  method  to  this  problem  by  Vaughn,  Oberkampf,  and  Wolfe. 
However,  no  attempt  will  be  made  to  be  complete  in  the  comparisons,  but  rather 
representative. 

For  0  <  Re  <  1.0,  f  =  0.1  and  A  =  3.0,  the  variation  of  C|_5|v|  with  Re  is 
linear  according  to  results  in  Figure  12;  the  maximum  discrepancy  in  C|_5|,^  is 
0.3  X  10-^  at  Re  =  1.0. 

In  Figure  13  C|_5f^  vs  Re  is  shown  for  f  =  0.1  and  A  =  3.0.  Over  the 

range  1  <  Re  <  100,  there  is  a  barely  discernable  maximum  at  Re  =  50, 

approximately.  For  Re  >  100,  Cl5|v|  decreases. 

The  variation  of  with  f  is  given  in  Figure  14  over  the  range  of 

0  <  f  <  1.1  for  the  spatial  eigenvalue  calculations  and  .05  <  f  <  0.9  for  the 
Strikwerda  results,  Kq  =  2°,  both  for  Re  =  10.0  and  A  =  3.0.  In  practice  f  < 
1.0  but  note  that  C^^f^  <  0  for  f  =  1.1.  The  two  sets  of  results  are  essen¬ 
tially  the  same,  which  tends  to  support  the  validity  of  each.  The  maximum 
difference  is  .0025  at  f  =  .6.  The  existence  of  a  maximum  in  this  curve  is  of 
some  interest.  When  the  first  measurements  of  C|_5f^,  from  yaw  growth  rate, 

were  made  at  "small"  Re,  a  linear  variation  of  Cl5|^  with  f  was  found,  but  only 

over  the  range  0  <  f  <  0.15.  It  was  speculated  that  for  "small"  Re  a 

Stewartson  resonance^  does  not  occur.  Figure  14  shows  that  resonance  does 
occur  at  Re  =  10,  albeit,  highly  damped.  Of  course,  Stewartson' s  theory  was 

formulated  for  an  inviscid  fluid  and  its  extensions  apply  to  "large"  Re  so 

that  quantitative  results  from  these  could  not  be  expected  to  apply  at  Re  = 
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A  comparison  of  calculated  side  moment  coefficient  vs  f  with  some  experi¬ 
mental  measurements  of  D'Amico  is  shown  in  Figure  15  for  Re  =  21.5  and 
A  =  1.042.  The  experimental  results  for  -  2°  were  deduced  from  yaw  growth 

rate  measurements  on  a  gyroscope;  according  to  Nusca  and  D'Amico, i®  the  Clsm 

and  f  are  in  S.  The  calculations  show  that  is  linear  for  0  <  f  <  .04  , 

approximately.  The  experimental  results  for  the  two  larger  values  of  f  lie  on 

the  calculated  Clsm  vs  f  line,  essentially.  The  reason  for  the  discrepancy  at 

t  =  0.045  is  not  known  although  it  is  pointed  out  in  Reference  19  that,  as 

Clsm  decreases,  the  error  in  its  measurement  increases;  no  error  bounds  were 

presented. 

A  similar  comparison  is  shown  in  Figure  16  for  Re  =  133  and  A  =  1.486. 
The  experimental  data  is  for  Kq  <  5°;  Clsm  deduced  from  yaw  growth  rate 

data.  The  calculations  show  that  Clsm  linear  for  0  <  f  <  .03,  approxi¬ 
mately.  The  gyroscope  data  are  uniformly  below  the  calculated  results;  no 
error  bounds  were  presented.  If  extrapolated,  the  trend  of  the  data  implies 
Clsm  =  O  at  f  =  0.03,  approximately,  and  negative  values  for  f  <  0.03.  In  all 

cases  computed  thus  far,  Clsm  positive  for  0  <  f  <  1. 

Next,  Clsm  vs  a  is  presentea  in  Figure  17  for  Re  =  10  and  f  =  0.1.  The 

variation  is  rapid  in  the  neighborhood  of  A  =  1  but  for  3  <  A  <  15,  Clsm 

creases  by  only  3%.  The  results  in  Figure  17  can  be  used  to  assess  the  assump¬ 
tion  of  infinite  aspect  ratio  in  Herbert's  theory. 3  For  Re  =  10  and  f  =  0.1, 
his  theory  gives  Clsm  '  which  is  quite  different  from  results  in  Figure 

17.  Since  the  contributions  of  Eqs.  (4.9a,b,d,)  to  Clsm  identically  zero 

according  to  his  theory,  another  point  of  view  is  to  compare  his  result  with 
that  obtained  from  only  Eq.  (4.9c),  using  the  present  results;  this  is  -0.184 
for  A  =  15,  Re  =  10,  and  f  =  0.1.  The  presence  of  the  endwalls  appears  to 
have  a  significant  effect  on  the  contribution  to  Clsm  shear  force  on 

the  sidewall,  even  for  A  =  15. 

The  variation  of  Clsm  ^  ’3  quite  different  for  Re  =  1000,  f  =  0.1  as 

shown  in  Figure  18.  The  calculations  using  LS  and  COL  both  produce  resonant 
response  curves  but  the  computed  values  of  Clsm  tl^an  in  thj 

comparisons  shown  in  Figures  4,  5,  and  6;  near  the  maximum  of  Clsm»  which 

occurs  at  A  =  1.12  for  COL,  they  differ  by  IIX.  For  Re  =  1000,  the  limiting 
values  of  Clsm  determined  as  accurately  as  for  smaller  Re  or  for  Cp; 

the  numerical  process  responsible  for  that  is  not  understood.  The  only 
guidance  as  to  which  of  the  calculations  is  the  more  accurate  is  provided  by 
the  third  set  of  results  in  Figure  18,  labeled  asymptotic,  computed  using  the 
theory  of  References  4  and  13  which  is  asymptotic  for  Re  +  ».  The  comparisons 
of  results  from  tnat  theory  with  measurements  of  yaw  growth  rate  given  in 
Reference  13  indicate  that  in  the  neighborhood  of  Re  =  1000,  f  =  0.1  and  A  =  1 


this  asymptotic  theory  is  reliable,  but  the  conclusion  is  not  definitive.  The 
differences  between  the  COL  and  asymptotic  results  are  less  than  the  corres¬ 
ponding  LS  differences  but  the  situation  is  too  tentative  to  draw  conclusions 
concerning  the  accuracy  of  either  method.  Both  LS  and  COL  results  should 
tend,  for  increasing  Re,  to  those  of  the  asymptotic  theory  for  decreasing 
Re.  Qualitatively,  this  tendency  is  illustrated  in  Figure  18  for  Re  =  1000. 

Despite  the  difference  just  discussed,  the  existence  of  the  resonant 
response,  wherein  the  cavity,  i.e. ,  liquid-filled  cylinder,  is  tuned  by  aspect 
ratio,  is  shown  in  Figure  18.  A  resonant  response  is  also  obtained  in  the  Cp 

vs  A  results;  the  maximum  Cp  is  at  A  =  1.12.  Usually,  resonant  response  tuned 

by  f  is  investigated  in  this  problem;  but  the  possibility  of  tuning  by  A  is 
not  surprising  considering  oscillatory  systems,  in  general.  For  fixed  A  = 
1.12  and  varying  f  calculations  by  the  asymptotic  method  give,  in  addition  to 
the  maximum  of  Cl^h^i  at  f  =  0.1,  a  secondary  peak  at  f  =  0.5. 

An  attempt  was  made  to  compare  the  moment  coefficients  calculated  by  this 
linear  theory  to  those  obtained  from  the  large  scale  numerical  metnods  in 
which  linearity  in  Kq  is  not  assumed.  In  addition  to  the  finite  difference 

method  of  Reference  6.  a  finite  element  method  was  recently  applied  to  this 
problem  by  Rosenblat^^  from  which  a  limited  number  of  results  is  available. 
Results  for  moment  coefficient,  with  Kq  =  20°,  are  presented  in  Table  1  from 

these  methods  and  the  spatial  eigenvalue  method.  The  test  was  frustrated 
because  the  results  using  the  methods  of  References  6  and  20  differ  signifi¬ 
cantly,  as  seen  in  Table  1.  Therefore,  the  effects  of  departure  from  linear¬ 
ity  at  Kq  =  20°  cannot  be  judged.  In  Table  1  Reynolds  number,  frequency  and 

moment  coefficient  are  given  in  the  aeroballistic  reference  frame.  S',  used  in 
References  6  and  20,  and  denoted  by  Re',  t',  and  0^5^';  the  transformations 

relating  these  to  the  corresponding  quantities  in  S  were  given  earlier.  The 
comparisons  are  made  in  S'  because  fewer  calculations  were  required. 

TABLE  1.  Comparison  of  C|_5|y|'  by  the  Spatial  Eigenvalue, 

Finite  Difference  and  Finite  Element  Methods; 
t'  =  0.1670,  A  =  4.29,  Kq  =  20°. 


Ref.  6,7 


Reference  20 


Spatial  E.V. 


5.90 

11.42 


0.03172 

0.4732 


0.03082 

0.04378 


0.03104 

0.04785 


The  differences  between  the  spatial  eigenvalue  results  and  those  of  Reference 


20  are  0.7%  and  10%  for  Re'  = 
spatial  eigenvalue  results  and 
differences  are  2.2%  and  1.1% 
between  these  methods  for  =  2 


5.90  and  11.42,  respectively.  Between  the 
those  from  Reference  7,  the  corresponding 
which  are  no  greater  than  the  differences 
’.  Although  a  wider  range  of  parameters  must 


be  investigated  to  arrive  at  a 
nonlinear  effect  is  small  up  to  Kq 


definitive  conclusion,  it  appears 

=  20°. 


It  has  been  demonstrated  that  the  linearized  Navier-Stokes  equations  with 
boundary  conditions  appropriate  to  a  spinning,  coning  cylinder  can  be  solved 
by  eigenfunction  expansions  over  a  wide  range  of  the  governing  nondimensional 
parameters;  Re,  f,  and  A.  By  expanding  the  complete  solution  in  powers  of  a 
small  parameter  related  to  coning  angle,  the  nonlinear  problem  can  also  be 
solved  by  eigenfunction  expansions.  The  limitations  of  the  linear  theory  could 
be  judged  by  comparison  with  finite  difference  or  finite  element  methods,  but 
this  is  not  yet  feasible. 

The  method  presented  here  is  in  the  class  of  spectral  methods.  Another 
spectral  method  using  Chebyshev  polynominals  has  been  developed  by  Herbert;2i 
results  from  his  method  are  not  yet  available  for  comparison.  Many  different 
solutions  for  the  fluid  motion  in  a  spinning,  coning  cylinder  are  now  availa¬ 
ble.  Judgments  on  the  relative  merits  of  these  should  be  made.  To  a  limited 
extent,  this  was  done  in  this  report. 

Certain  deficiencies  in  the  spatial  eigenvalue  method,  for  Re  >  500,  were 
noted.  Ideas  for  removing  these  are  being  exploited. 
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Figure  2.  The  first  15  eigenvalues  for  the  Stokes  limit  and  for 
Rg  =  100.  f  =  0.1. 
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Fiqjre  5.  Pressure  coefficient  at  r  =  0.667  vs  number  of 
eigenvalues  for  Re  =  100,  f  =  0.1,  and  A  =  3.0 
calculated  by  LS  and  COL. 
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Figure  6.  Pressure  coefficient  at  r  =  0.667  vs  number  of 
eigenvalues  for  Re  =  1000,  f  =  0.1,  and  A  =  3. 
calculated  by  LS  and  COL.  _ 
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Figure  7.  Comparison  of  Cp(r,A)  from  the  Strikwerda,  - 

and  spatial  eigenvalue  methods  for  Re  =  1.0, 
f  ^  0.1,  and  A  ^  3.0. _ 
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Pressure  coefficient,  at  x  =  A  and  two  values  of  r,  v 
for  Re  =  10  and  A  =  3.0,  according  to  the  spatial  eig 
value  method  and  the  asymptotic  approximation  (4.6). 
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Figure  12. 


Moment  coefficient  vs  Re  for  f  =  0.1,  A  =  3.0 

according  to  the  spatial  eigenvalue  method  and 

Strikwerda's  method,  K  =  2°,  0  <  Re  <  1.0. 
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Figure  15.  Moment  coefficient  vs  f  for  Re  =  21.5,  A  =  1.0 
according  to  the  spatial  eigenvalue  method  and 
measurements  from  Ref.  19  for  K.  =  2°. 


Moment  coefficient  vs  f  for  Re  =  133,  A  =  1.486  according  to  the  spatial 
eigenvalue  method  and  measurements  from  Reference  19  for  K  =  2°. 
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Figure  18.  Moment  coefficient  vs  A  for  Re  =  1000,  f  =  0.1,  0.9  5  A  <  1.3 

according  to  the  spatial  eigenvalue  method  using  LS  and  COL  and 
including  results  according  to  the  asymptotic  method  of  Ref. 
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